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Abstract 

We generalize the concept of nonlinear periodic structures to systems 
that show arbitrary spacetime variations of the refractive index. Non- 
linear pulse propagation through these spatiotemporal photonic crystals 
can be described, for shallow nonstationary gratings, by coupled mode 
equations which are a generalization of the traditional equations used for 
stationary photonic crystals. Novel gap soliton solutions are found by 
solving a modified massive Thirring model. They represent the missing 
link between the gap solitons in static photonic crystals and resonance 
solitons found in dynamic gratings. 

The ability to manipulate spectrally and temporally optical pulses has been 
a longstanding goal in modern science and technology, and periodic media have 
the potential of engineering light propagation to an unprecedented degree pQ. 
Systems possessing a periodic modulation of refractive index, in which pho- 
tonic bandgaps (PBGs) form at or near a multiple of the Bragg frequency (or 
wavenumber) , have been used in a wide range of applications, including disper- 
sion compensators and optical filters pQ. When Kerr nonlinearity is considered, 
new effects come into play, such as optical bistability, pulse compression, optical 
switching and soliton formation pQ. 

Stationary gap solitons (GSs) living in the frequency bandgap of a ID pe- 
riodic medium were first investigated in 1987-1989 in a series of fundamental 
papers [H 02 |3|, and found experimentally in 1996 [5]. Solitons living in a 
wavenumber bandgap of the so-called dynamic gratings (i.e. a traveling-wave 
periodic index change) were also investigated in Refs. [9l[TTJ[l0] by using coprop- 
agating beams, where the complementarity between the two kinds of bandgap 
was evident. Since then, there has been an exponential increase in the number 
of studies and practical applications of GSs. An intriguing possibility is the 
storage of optical pulses in the form of zero velocity GSs followed by release 
from the structure at a controllable delay pQ. 

Much less attention, however, has been devoted to the physics of nonstation- 
ary periodic media, such as dielectric structures showing temporal variations of 
the refractive index [7], which have the potential to dramatically enhance 
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the degree of spectral control over light pulses by periodic media thanks to the 
new temporal degree of freedom [BJ. In Ref. [7J we derived the transfer matrix 
T for plane waves scattered by the sharp boundary associated with a medium 
with time-varying refractive index, which must be distinguished from a mov- 
ing interface in that the medium itself is immobile [8] . The knowledge of T 
for a single boundary allowed us to construct a theory for more complicated 
nonstationary dielectric objects. In particular in [7] we introduced the impor- 
tant concept of spatiotemporal photonic crystal (STPC), which is a grating that 
shows a well-defined periodicity of the refractive index along a certain direction 
of the spacetime plane (z, ct) (z is the longitudinal spatial coordinate, t is time 
and c is the speed of light). This periodicity gives rise to PBGs in a mixed 
frequency-wavenumber space, the mixing being regulated by an angular param- 
eter 9, which we shall see it is related to the apparent velocity of the layers in 
the spacetime plane. 

In this Letter we extend the linear theory formulated in Ref. [7J to non- 
stationary gratings with Kerr nonlincarity, demonstrating the existence of self- 
localized solutions in the mixed frequency-wavenumber bandgaps of STPCs, 
thus showing that the conventional concept of GS (see Refs. [3J [U ® ITU] ) 
must be extended to encompass general spacetime variations of the refractive 
index. 

Let us consider an electromagnetic wave, with its electric and magnetic fields 
E = (E(z, t), 0, 0) and B = (0, B(z, t), 0) linearly polarized along the x and y 
directions respectively. E and B depend on z and t only, because we assume 
conditions of normal incidence, so that any change of the time-dependent re- 
fractive index occurs along z. The linear polarization of the medium is given 
by Pi, = xl(z, i)E = [n(z, t) 2 — 1]E, where \l is the linear susceptibility of the 
(non-magnetic) medium and n(z,t) is the linear refractive index, which is as- 
sumed for simplicity to be real and frequency independent, and possessing for the 
moment an arbitrary dependence on z and t. Maxwell's equations for E and B 
are dt{eE)+cd z B+d t P NL = 0, d t B+cd z E = 0, where e(z, t) = = n(z,t) 2 , 

and Pnl is the Kerr nonlinear polarization, Pnl = XnlE 3 , with xnl constant. 
Here and in the following we use the Heaviside-Lorentz units system, see Ref. 
0. 

By deriving the first of Maxwell's equations with respect to t, and using the 
second one to eliminate B, we obtain the nonlinear wave equation for a space- 
and time-varying refractive index: 

ed 2 E - c 2 d 2 z E + {8 2 e)E + 2(d t e)(d t E) + d 2 P NL = 0. (1) 

It is now convenient to introduce two new variables, rotated by an angle 9 £ 
[— 7r/2,7r/2] in the (z, ct) plane: p = cos(9)z — sin(9)ct,q — sin(6>)z + cos(9)ct. 
This spacetime rotation is analogous to the Lorentz transformations in special 
relativity, with the essential conceptual difference that in our case the associated 
dimensionless velocity (3 = tan(#) can assume values in the range < |/3| < oo, 
and it is not limited by c, see Ref. [7] and references therein, p will be chosen 
to correspond to the direction parallel to the periodicity of the STPC, while q 
will be orthogonal to p. 9 > implies that the boundaries of the STPC are 
moving towards light, while for 9 < the boundaries are moving away from the 
incident pulse. A generalized plane wave propagating in the (p, q) space has the 
form Vl 1 = ^Qexp(ikp — iujq), where is a constant amplitude, and k and Q 
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Figure 1: (a) Space- like, conventional static photonic crystal (8 — 0), for which 
p — z and q = ct. and are respectively refractive indices and widths of 
the two types of layers, (b) Time-like photonic crystal, 9 = ±7r/2, for which the 
refractive index changes periodically in time only (p = —ct), and q = z. ct\^ 
are the durations of the layers, (c) Intermediate case < \9\ < ir/2. A p is the 
period of the structure along the grating direction. Axis ct and z are indicated 
in a circle. 



are the wavenumbers associated to the p and q directions respectively, k and u> 
are linked by the rotation k = cos(8)k + sm(9)oj/c and Hi — cos(9)oj/c — sin(9)k, 
where k and co/c are the wavenumbers associated to the original physical plane 
(z,ct). Figure [T] shows the geometrical meaning of axes p and q for three rep- 
resentative STPCs of fundamental importance (see caption). It is evident from 
the above definitions that the case 9^0 corresponds to layers arranged peri- 
odically along z (p — > z), and the plane wave derealization direction lies along 
ct (q — > ct). This corresponds to the traditional time-independent photonic 
crystal, see Fig. [IJa). Being the crystal invariant with respect to translations 
along ct, we name this structure a space-like STPC. More interesting is the sec- 
ond limiting case, when 9 —> ±ir/2, shown in Fig. |Tfb) . From the definitions 
of p and q, we have p — > ^fct and q — > ±z, so that plane waves will be delo- 
calized along z, and all the variations of refractive index occur in time only. 
In analogy with the previous nomenclature, we name this structure a time-like 
STPC, an example of which is the dynamic grating [TTJIH]. The intermediate 
cases when < \9\ < tt/2, which are the main focus of this Letter, are displayed 
schematically in Fig. [He). 

In principle, integration of Eq.([T]) is all one needs to completely solve the 
problem of nonlinear pulse propagation in any kind of nonstationary dispersion- 
less structure. However, one can gain important analytical insight by consid- 
ering a cosinusoidal shallow grating described by the dielectric function e(p) = 
n\ + fi[exp(ikBp) + c.c], where Uq is the square of the average linear refractive 
index, and /x <C Uq. Here, ks represents the equivalent of the Bragg wavenum- 
ber along p. The Bragg condition, which is the phase-matching condition be- 
tween the optical and the grating wavenumbers, is given by ks — k + — k~ . 
Due to the fact that only two Fourier modes are present in the above ex- 
pression for e(p), we can assume that only two optical modes strongly con- 
tribute to the propagation dynamics, which leads us to the expansion E(p, q) = 
[F(p,q)e lk+p ~ lC:}q + B(p,q)e ik ~P- lQ "i + c.c.]/2, where F and B are envelopes of 
respectively the forward and backward components of the electric field, k^ are 
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the wavenumbers along p for F and B respectively, and u> is the wavenumber 
along q, which is common for both components. Note that: (i) here the terms 
'forward' and 'backward' are not in general associated to the spatial motion 
of the modes (unless 6 = 0), but rather to the more general motion along the 
p-direction, and (ii) the F and B modes will generally have different linear 
wavenumbers \k \ along p. The dispersion relation between k and Cj is (see 
also [7]) = [sin(0) ± no cos(#)]u)/[cos(#) =p n osin(#)], which is not valid for 
= ± arctan(l/no). For those angles, either k + or fc _ diverge and the wavenum- 
ber along the p-direction is not defined. Physically this is due to the fact that 
for 6 < — arctan(l/no) the layers boundaries of the STPC are changing faster 
than the speed of light in the medium (c/no). 

After substituting the expansion for E(p,q) into Eq.fl]), a slowly- varying 
amplitude approximation (SVEA) in p and q is performed: \d 2 ip\ <C |fc ± 9 P V| ^ 
[(k^) 2 ^, \d 2 ip\ <C \Qd q ip\ <§; \oj 2 ^\, where ip is either F or B, and similar rela- 
tions are valid for the terms containing the mixed derivative d p d q . The following 
two spatiotemporal coupled mode equations (STCMEs) are obtained: 

+ n cos(0) + sH0) + n + T + = 

cos(w) — n sm(0) [cos(0) — no sm(0)\ z [cos(#) - n sm(0)\ z 

n cos(g)-sin(g) . . " . — F + . j- T , , mi2 (2|F| 2 + |B| 2 )i? = 0, (3) 

p cos(0) + n sm(0) 9 [cos(6») + n sm(6»)] 2 [cos(6>) + n sm(6»)] 2 v 11 1 ; w 

where T = (3wx./vl)/(871o) is the nonlinear coefficient, and k = w/i/(2no) is 
the grating coupling constant. Again, the equations have a singular character 
for = ± arctan(l/rio). Eqs.([2][3|) represent the first central result of this Let- 
ter. Let us now perform the following scaling, which is well-defined for ^ 
arctan(l/no) and ^ — arctan(no): r = q/qo, £ = p/po, f = F/Aq, b = B/Aq, 
withpo = K~ 1 [cos(9)—n o sm(0)] 2 , q a = [n cos(9) + sin(9)]p /[cos(9) — n sin(6>)], 
Fo = (re/r) 1 / 2 . With this, equations (j^EJ) are reduced to the following two di- 
mensionless equations: 

t(d i +d T )f + b+(2\b\ 2 + \f\ 2 )f = 0, (4) 
i (-d e + pid r ) b + p 2 f + p 2 (2|/| 2 + \b\ 2 ) b = 0, (5) 

where pi(9) = [no cos(0) — sin(#)] [cos(#)— no sin(#)]/{[cos(6>)+no sm(9)} [no cos(#)+ 
sin(0)]} and p 2 {6) = {[cos(0) - n sin(0)]/[cos(0) + n Q sin(<9)]} 2 . In Figure^a) 
coefficients are shown as a function of 9 for an average refractive index 
no = 3. Note that although p 2 is always positive in the range 9 e [— 7r/2, 7t/2], 
Pi becomes negative for arctan(l/no) < [9[ < arctan(no), and shows two di- 
vergences for negative angles at 6 = — arctan(no) an d at 9 = — arctan(l/no). 
Also pi = p2 = 1 for the limiting cases 6 = and 9 = tt/2, but in general these 
parameters can strongly differ from unity. 

Let us now discuss the most important linear property of Eqs.(j4][5]), namely 
the PBG in the 0-rotated frequency-wavenumber space. Substituting {/, 6} = 
^f,b exp(ifc'£ — iu'r) into Eqs. (|3EJ) , and neglecting the nonlinear terms, we 
readily obtain uj' 12 (0) = {( Pl -l)k' ±[(l+ Pl ) 2 k' 2 +Ap 1 p 2 } 1 / 2 } / (2 Pl ), after which 
we perform an inverse rotation back to the original dimensionless frequency- 
wavenumber space, i.e. k" = cos(9)k' — sm(0)ui' , uj" = sm(0)k' + cos(0)ui' . In 
Figure [^d,e,f) the bandstructure Lo" 2 {k") for three different cases (0 = 0, = 
1.3 and = 7r/2) is plotted, explicitly showing the passage from the frequency 
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Figure 2: (Color online) (a) ^-dependence of coefficients pi.2 for the parameter 
range 9 6 [— arctan(l/no),7r/2], and for 6 G [—n/2,0] (inset), (b) Contour 
plot of Itot = |/| 2 + \b\ 2 for numerical propagation of an initial STGP (input 
profile shown inset, blue line is |/| and red line is with parameters v = 0, 
9 = 1.3 and 6 = tt/2, without grating (k = 0). Propagation length is r = 7. 
(c) Same as (b) but with grating coupling, (d) Bandstructure in space (k",u>") 
as calculated by using (JHEJ), for 9 = 0, (e) for 9 = 1.3 and (f) for 9 = tt/2. In 
(c,d,e), light gray regions indicate the w"-bandgap, and the dark gray regions 
the fc"-bandgap. In all cases uq = 3. 



bandgap [9 — 0, Fig. [H(d)] to the wavenumber bandgap [8 = tt/2, Fig. [H(f)], 
passing through a region in which the two kinds of bandgap coexist [9 = 1.3, 
Fig. He)]. 

We now proceed to analyze the symmetries and the GS solutions of Eqs.((3]- 
[H]). One can derive Eqs.(|3][H]) from the following Hamiltonian density: H = 
(bf* + fb*) + 2\f\ 2 \b\ 2 + (|/| 4 + |6| 4 )/2 - M f + M b / P 2, where the star indicates 
complex conjugation, and where M.c, = *(C^?C* — c.c.)/2 = Im{C*9^C} is the 
momentum density of the generic field £ = {/, b}. The dynamical equations are 
written as iJAld T g + 5H/Sg^ = 0, where M = diag(l, p\/ P2), J = diag(l, — 1) 
is the symplectic matrix, g = [/, b] is the field vector, and dagger indicates 
hcrmitian conjugation. Moreover, the variational derivative is given by 5/(S() = 
d/(d() — d^[d/d(d^()}, see also Ref. [H]. We anticipate that £ will correspond 
to the localization coordinate of the soliton solutions, and r to the evolution 
coordinate. With this in mind, one can find the total Hamiltonian by integrating 
H over f , H = Hd£ = [H]±™. H is an integral of motion, i.e. d T H = 0. H 
does not depend on the variable ^ explicitly, leading to the conservation of total 
momentum: M to t = [Mf + (pi/ P2)-Mb]-^, i-e. d r M to t = 0. H is also invariant 
with respect to the 'gauge transformation' / — > / exp(i</>), b — > 6exp(i0), leading 
to the conservation of the quantity P = [|/| 2 + (pi/p2)|^| 2 ]l^, and d T P = 0. 
The number of integrals of motion of the dynamical system determined by (f5][5]) 
(with the exclusion of H) is closely related to the number of internal parameters 
of the corresponding soliton families |13j . Therefore the family of localized 
solutions living inside the li-bandgap of a STPC are represented by two internal 
parameters. This is well-known for solitons living in the frequency bandgap of 
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a static photonic crystal (9 — 0) [3J S] , and for GSs living in the wavenumber 
bandgap (9 = ±tt/2) [HUH], but our analysis extends this result for arbitrary 
values of 9. 

In order to find analytical localized solutions of Eqs.fjUE]), let us now consider 
the following different set of coupled equations for two new fields ipf and tpb- 

i(d e + d T )ipf+iJj b + \iJj b \ 2 iJjf = 0, (6) 
i (-<% + pid T )ip b + p 2 ip f + p 2 \ip f \ 2 if) b = 0. (7) 

In analogy with the extensively studied Massive Thirring Model (MTM) [HI H] , 
we name Eqs.piS]) the modified Massive Thirring Model (mMTM). The MTM 
is a particular case of the mMTM, with p\ = p 2 = 1, and it is known to 
be integrable [T3]. The mMTM solitons will automatically provide analytical 
soliton solutions to the original equations Eqs.(@][5]). Let us operate the Galilcian 
shift r' = t — V£, C = £, and the rescalings f = t'/tq and £ = By choosing 
tq = 1 + V and V = (1 — px)/2 we can therefore scale p\ away from Eqs.|6|[7]). 
It is now possible to find the following analytical soliton solution: 

</>/(£» = (l/A)Vo, MLt) = -AiIj , (8) 

1/2 

where ipo = sin((S)sech (© — iS/2) exp(zcr), 7 = [^2/(1 — « 2 )] , A = [^2(1 — 
v)/(l + v)] 1 ^, 9 =_7sin(*)(f- ut) = 7 sin((5){[l + uV/(l - V)]£ - Vt/(1 - V)} 
and a = 7 cos(5)« - f) = 7 cos(5){[u + V/(l - V)]£-t/{1-V)}. < (5 < vr is 
a parameter (also called the soliton charge) which measures the detuning from 
the bandgap center (S — ir/2), and — 1 < v < +1 is the soliton relative velocity. 

We now attempt to express the general localized solutions of Eqs.(|3][5]) in 
terms of mMTM solitons, Eqs.©, by using the ansatz: {/, b} = aijjffi{£,, t) exp[w7(0(£, t))]. 
Substituting into Eqs.(j4][5|) and using ([8|), we obtain two equations for rj' = 
Dqt]. The consistency condition between them determines the value of a: 
a = {[2(1 - v 2 )p 2 ]/[(l + v) 2 + 4(1 - v 2 )p 2 + (1 - vfp 2 ]} 1 ' 2 , which in turn 
is used to solve the ODE for 77(8), obtaining the solution: e"'( e ' = {[1 + 

e i5+2Sj n e ii + e 2e^[(l+v) 2 ~(l~v) 2 p 2 2 ]/[(l+v) 2 +4(l-v 2 )p 2 + (l-v) 2 p 2 2 }^ Thig completes 

the information necessary to find the two-parameter family of localized solu- 
tions, i.e. the spatiotemporal gap solitons (STGSs), for the STCMEs given by 
Eqs.(@][5]), which represents the second central result of this Letter. The inten- 
sity ratio r between / and b is given by r = \f\ 2 /\b\ 2 = (1 + v)/[(l — v)p 2 ], so 
that / and b for the zero-velocity solitons (v — 0) do not have in general equal 
amplitudes [r(v = 0) = l/p 2 ]- Figure [3] shows contour plots of the soliton total 
intensity I to t — \f\ 2 + \b\ 2 when changing 9 [Fig. [3(a)], S [Fig. [2(b)] and finally 
v [Fig. He)]. 

Our analytical solution is also confirmed by direct numerical integration of 
Eqs.([3][5]), performed using a split-step Fourier method with a 4th order Runge- 
Kutta algorithm. Figure [2](b,c) shows the propagation of a STGS with parame- 
ters v = 0, 6 = 1.3 and 6 = ir/2, which lives in the center of the mixed bandgap 
displayed in Fig. Ute), for a propagation of t = 7. Fig. [U^b) shows that when 
the grating is absent (k = 0) the two components separate and do not inter- 
act, while Fig. Oc) shows the undisturbed soliton propagation at zero relative 
velocity in presence of the spatiotemporal grating. Surprisingly, it is seen from 
the numerical simulations that quasi-adiabatic variations of 9 during propaga- 
tion, which thus change dynamically the background spatiotemporal grating, do 
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Figure 3: (Color online) Contour plots of the total intensity profile I to t = 
|/| 2 + |6| 2 for the STGS: (a) as a function of 9, for v = and 5 = ir/2 [note the 
divergence located at 8 — arctan(l/no), and the vanishing intensity in corre- 
spondence of 6 = — arctan(l/no)], (b) as a function of 6, for 8 — 1.3 and v = 0, 
and (c) as a function of v, for 8 = 1.3 and 6 = ir/2. In all cases no = 3. 



not destroy the STGS, due to prompt pulse reshaping. This structural stability 
makes STGSs very attractive for storing, slowing down, converting and releas- 
ing optical energy in a controlled way, which may have profound implications 
for optical communications and quantum information processing [15] . 

In conclusion, in this Letter we have derived a set of CMEs that allow to 
describe nonlinear pulse propagation in a shallow grating with space-time vari- 
ations of the refractive index. This structure generally possesses a bandgap in a 
rotated frequency-wavenumber space, where new GSs have been found analyti- 
cally by solving an associated mMTM. Our formulation considerably generalizes 
the current theoretical understanding of periodic media to time-dependent re- 
fractive index. Future works will include the bifurcation and stability analysis 
of STGSs, and the natural extension of the theory to dispersive media. 

We acknowledge financial support from the UK Engineering and Physical 
Sciences Research Council (EPSRC), Science Foundation Ireland (SFI) and the 
Irish Research Council for Science, Engineering and Technology (IRCSET). 
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